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Abstract — We describe a novel approach to morphology and 
anisotropy analysis of complex spatial structure using tensor- 
valued Minkowski functionals, the so-called Minkowski tensors. 
Minkowski tensors are generalizations of the well-known scalar 
Minkowski functionals and are explicitly sensitive to anisotropic 
aspects of morphology, relevant for example for elastic moduli 
or permeability of porous materials. We provide explicit linear- 
time algorithms to compute these measures for three-dimensional 
shapes given by triangulations of their bounding surface, includ- 
ing triangulations obtained from experimental gray-scale image 
data by isosurface extraction. Eigenvalue ratios of Minkowski 
tensors provide a robust and versatile definition of intrinsic 
anisotropy. This analysis is applied to bio-polymer networks 
under shear and to cellular complexes of dense bead packs. 

Index Terms — Mathematical morphology, anisotropy charac- 
terization, porous and cellular structures, Minkowski functional 

The morphology of complex spatial microstructures is 
often classified qualitatively into types such as cellular, 
porous, network-like, fibrous, percolating, periodic, lamellar, 
hexagonal, disordered, fractal, etc. Various quantitative 
measures of morphology have been defined often applicable 
to one specific type only, for example moments of the 
distributions of angles of tangent vectors with a fixed 
specified direction as anisotropy characterization of a network 
structure. Apart from the concept of correlation functions, 
few measures are defined sensibly and robustly for all types. 
In this article, we describe the class of Minkowski tensors 
(MT) that apply generically to any type of structure, which 
contains two or more phases separated by a well-defined 
interface (e.g. porous media, foams, trabecular bone, granular 
material in Fig. 1). The MT are defined as integrals of powers 
of normal and position vectors and surface curvatures (or 
curvature measures). Because of their tensorial nature they 
are explicitly sensitive to anisotropic and orientational aspects 
of spatial structure. Robust indices of intrinsic anisotropy 
and alignment can be derived. Figure 1 shows examples 
of systems where subtle anisotropy of the spatial structure 
influences the physical properties and to which the analysis 
of this article is applicable. The scope of this article is to 
supply a thorough description of an algorithm to compute MT 
of bi-phasic structures, when the phase-separating interface 
is given by a triangulation. Two examples illustrate the 
Minkowski tensor analysis. For further applications of MT 
for pattern analysis we refer to refs. [l]-[5], [5]-[7]. 

Scalar Minkowski functionals (MF) are well established 
as succinct descriptors of morphology and spatial structure 



for various physical processes [10]. These integral geometric 
measures have been applied to disordered porous materials 
[11] and are relevant to flow phenomena therein [12], to nano- 
scale microstructures in copolymers [13], to the dewetting 
dynamics of thin films [14], and to Turing patterns [15]. They 
have also been shown to be the most pertinent morphological 
quantities on which the thermodynamic properties of simple 
fluids near curved solid interfaces depend [16], [17]. The 
mathematical theory of MF and their generalizations has 
been comprehensively developed in the context of integral 
geometry [18]— [21], with several aspects shared also with the 
discipline of mathematical morphology [22]-[24]. 

MF as scalar quantities are not sensitive to features of 
the morphology which relate to orientation or directional 




(a) Copolymer Film in Field (b) Metal Foam 




(c) Trabecular Bone (d) Granular Material 

Fig. 1. Examples of systems with anisotropic spatial structure, (a) A 
microphase-separated copolymer film aligns under the influence of an external 
electric field (adapted from [8], reproduced by permission of The Royal Soci- 
ety of Chemistry), (b) Closed-cell metal foam (image courtesy of M. Saadatfar 
[9]). (c) Structure of trabecular bone (image courtesy Alan Boyde). (d) Packing 
of ellipsoids as systems for anisotropic granular matter. 
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anisotropy, since motion-invariance is one of their defining 
properties. Therefore they do not provide quantification of 
anisotropy that is relevant to study the direction-dependence of 
physical processes, such as elastic properties or permeability 
of anisotropic porous or microstructured materials or systems 
with external fields. This motivates their generalization to 
tensorial quantities. MT have already been shown to be the 
relevant morphological descriptors for a density functional 
theory of fluids of non-spherical particles [25] and of DNA 
conformations [26], and of a simple model for transport 
with molecular motors [27]. They have also been used, in 
2D, as morphology descriptors of arrangements of neuronal 
cells [28], galaxies [29], and Turing patterns [6]. The 
mathematical discipline of integral geometry has proven 
statements regarding continuity and completeness similar to 
the Hadwiger theorem in the scalar case [30]-[33]. However, 
an algorithm for the computation of the MT applicable 
to experimental 3D data - a prerequisite for their use as 
shape indices for experimental data - has thus far been lacking. 

A primary application of rank-2 MT is the quantitative 
analysis of the degree of intrinsic anisotropy of materials 
with complex spatial structure. Scalar measures of anisotropy 
are easily derived as eigenvalue ratios of the MT. Alternative 
methods for the characterization of anisotropy and alignment 
exist. Fourier transforms are a common way to characterize 
anisotropy, and have been applied e.g. for trabecular bone 
[34], for electrodeposited patterns [35], for fiber systems [36], 
and for structured polyethylene mats [37]. Related methods 
based on correlation functions are also known [38], [39]. 
Fourier methods that analyze the amplitude of the Fourier 
transform of a gray-scale image in polar coordinates can 
also quantify alignment, e.g. of copolymer films in electric 
fields [8]. Anisotropy indices derived from the normal vector 
distribution of a given shape, similar to the MT, have been 
used to describe the shape anisotropy of simulated 3D foam 
cells [40], [41] and liquid interfaces [3], [42]. An anisotropy 
measure applicable to porous media is derived from the 
directional variations of average chord lengths. For a binary 
composite, i.e. consisting of a solid and a void phase, a 
chord is a segment of an infinite straight line that is fully 
contained in one of the two phases. Analyses of chord lengths 
and the derived mean intercept length ellipsoid are used for 
the investigation of the microstructure of bone [43]-[48], 
see also ref. [49] for a comparison of anisotropy measures 
based on mean-intercept length, star-volume and star-length 
distributions. Deformations of cellular or granular material 
have recently been quantified using the so-called texture tensor 
C, defined as the sum Cy := Yl ^ih over a subset of link 
vectors 1 in the structure [50], [51]. The texture tensor can be 
used to characterize anisotropy, e.g. for Antarctic ice crystals 
[52] and liquid foam cells [53]. Further anisotropy measures 
are based on the Steiner compact [54], wavelet analysis [55], 
the orientation of volumes [56] or star- volumes [57]. Two- 
dimensional equivalents of the anisotropy measures discussed 
in this article have previously been used for the analysis of 
the shape of neuronal cells [28] and galaxies [29], and are 
discussed in detail in [6]. (Parts of the analysis of this article 




Fig. 2. (a) A body K with bounding surface dK; (b) a convex polytope P 
and its parallel body (or dilation) P tfcl B e ; (c) a subset of a topologically more 
complex body based on Craig Marlow's painting "White Spirits" [59]. The 
latter demonstrates a common experimental situation, namely that the given 
body K represents a finite subset of a larger body K + , here the body K 
clipped to the window of observation T, K = K + n T. K is assumed to 
be a representative subset of the larger body, which allows for the estimation 
of intrinsic shape features of K + . If only K, but not K + , is available for 
analysis particular care must be taken w.r.t. those bounding surface patches 
of K that result purely from taking the subset, i.e. those that are on the 
boundaries of the window of observation. 

represent the thesis work in ref. [58]). 

The paper is organized as follows: Section I provides 
an overview of the theory of MF and MT, based on their 
definition by surface integrals which is widely used in the 
physical sciences. To bridge the gap in notation between 
the physics and the mathematics literature this section also 
includes a discussion of the alternative definition based on 
measure theory. Section II describes algorithms to compute 
MT for triangulated bodies. Section III describes anisotropy 
measures derived from the MT and illustrates their application 
to two experimental data sets. The appendix provides analytic 
expressions for the MT for some simple geometric shapes. 

I. Definition and fundamental properties 

Scalar MF and MT can be used as shape measures to 
quantify the shape or form of an object. They can be defined 
in two largely equivalent ways. In the physical sciences, a 
definition based on curvature-weighted integrals of position 
or surface normal vectors over the object's bounding surface 
is commonly used, and forms the basis of the numerical algo- 
rithms derived in this article. An alternative, more fundamental 
definition is provided by the measure theoretic approach of 
integral geometry (see I-B). 

The object, also referred to as body, whose shape can be 
characterized by MT or MF is denoted by K. Assuming 
that K is a compact set with nonempty interior embedded in 
Euclidean space R 3 and is bounded by a sufficiently smooth 
surface dK, we define MF of K as 

W (K)=JdV and W„{K) = ^Jg v &A (1) 

K dK 

in space-dimension d = 3 and with v = 1,2,3. The scalar 
functions G v are G\ = 1, the mean curvature G2 = (fci + 
fc2)/2 and the Gaussian curvature G3 = k\ ■ &2 of the bounding 
surface dK (k\, ki are the principal curvatures on dK as de- 
fined in differential geometry); dV^ is the infinitesimal volume 
and dA the scalar infinitesimal area element. This definition 
naturally applies to both convex and non-convex bodies of 
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arbitrary topology with a sufficiently smooth bounding surface. 
The prefactor is chosen such that for a sphere Br with radius 
R the scalar MF are W V {B R ) = k 3 R 3 -" where k 3 = 4?r/3 
is the volume of the 3-dimensional unit sphere l . 

The MF W2 and W3 (in 3D space) are not properly defined 
by eq. (1) for bodies with sharp edges or corners, due to 
singularities of the mean and Gaussian curvatures G2 and 
G3. However, for convex bodies, consideration of a parallel 
or dilated body in the limit of vanishing thickness provides a 
robust definition, by use of the Steiner formula. The Steiner 
formula states that, for a given convex body K, the MF of the 
parallel or dilated body (K W B t ) are a polynomial in e with 
coefficients proportional to the MF of K; K e := K ttJ B e := 
{xi + X2 : xi G K, X2 G B e } is the parallel or dilated 
body of K (see Fig. 2). Specifically for the volume one finds 
W (KVB e ) = Wo(K) + 3Wi(K)e + 3W 2 (K)e 2 + W 3 (K)e :i , 
and more generally 



The MT of rank two are defined as 



W v {K\HB t ) = 



fj, — V x ' 



f-L — V 



(2) 



Sharp edges and vertices of K correspond to cylindrical or 
spherical segments on the bounding surface d{K y5 e ) of the 
parallel or dilated body. The bounding surface is sufficiently 
smooth for e > and the limit e ->■ of W V (K l±l B e ) 
converges to W V {K). It is further necessary to define MF and 
MT for certain non-convex bodies (with or without positive 
reach, see e.g. [21, Note 1 to Section 5.3 and refs. therein]); 
this is achieved below by exploiting an additivity relationship. 
A further discussion for non-smooth bodies can be found in 
section II. 

MT are symmetric tensors (that is, invariant under index per- 
mutation), which are generated by symmetric tensor products 
of position vectors x and normal vectors n of dK. The dyadic 
(or tensor) product of two vectors a and b is (a®b)ij := suhj. 
Let now a and b be symmetric tensors of rank r and s resp., 
then the symmetric tensor product is defined as 

(a©b) 4i ... Ws := 

T \| X] a <r(il) ' ' • a o-(^)b<T(i r+1 ) ■ • • b CT ( i? . + 3 ), (3) 



(r + s)\ 



where S r+S is the permutation group of r + s elements. We 
use for two tensors a and b the shorter notation a 2 := a© a = 
a <g> a and ab := a b. For example, if a and b are both 
vectors the symmetric tensor product is ij of the rank 2 tensor 
ab. 

'Other normalizations are common in the literature too. The kinematic 
formulae [20], [21] have particularly simple coefficients if the normalization 
M V (K) = K d _ l/ W v (K)/(K v K d ) is used, with the volume of the n- 



dimensional unit ball k„ 



7r"/ 2 /r(n/2 + 1). In the mathematical 
literature, in d-dimensional Euclidean space the normalization Vd- V {K) = 
W V (K)/ k v is frequently used, and the Vd- V are called the intrinsic 
volumes of K. In three dimensions, the set of MF thus consists of the volume 
Wo = Mo = V3, the surface area 3W\ = 8M1 = 2V2, the integrated 
mean curvature SW2 = 2tv 2 M2 = nVi, and the Euler characteristic 



W*'°(K) := Jx 2 dV 



K 
1 



(4) 



W(K) := -jG^ r n s dA. (5) 

dK 

with v = 1,2,3 and (r, s) = (2, 0), (1, 1) or (0, 2). For ease 
of notation, we set W^ s := for s > and W^ s := if 
v < or v > 3. For a three-dimensional body, this definition 
yields 10 MT (not counting the ones that vanish by definition 
for all bodies). 

MT of rank one (called Minkowski vectors) are defined by 
Wq'° := J K x dV and by := ± J gK x(L4 for v = 1, 2, 3. 
The prefactors are chosen such that, for a sphere centered 
at C, the so-called curvature centroids W^,'°/W v are equal 
to C. Note specifically that W^ 1 /Wq is the center of mass 
(assuming a solidly filled body of constant density). Formally, 
vectors W®- 1 proportional to J dR ndA for v = 1,2,3 are 
also defined, however they vanish for any body (with a closed 
bounding surface) [32]. 

MT are isometry covariant, that is their behavior under 
translations and rotations is given by 

w r ^{K\at) = J2( r \"wr p ' s (K) (6) 

W:' S (UK) - U r+s :W r v > s , (7) 

where K l±l t is the translation of K by the vector t, UK 
is a rotated copy of K, and U r+S denotes the corresponding 
rotation tensor for a rank-(r + s) tensor: 



U r+S : Wl 

ll,...,J r + s 



E 

jl,---,jr+s 



In this expression, Uij is the conventional orthogonal 3x3 
transformation matrix associated with the rotation U. 

For r = 0, equation (6) gives W^ S (K) = W^ S (K W t). A 
tensor that fulfills this relation for all K is called translation 
invariant. Genuinely translation covariant tensors fulfill (6) 
but not the translation invariance condition. For the sake of 
brevity, we will use the term translation covariant to denote 
specifically the genuinely translation covariant tensors. All 
W®' s are translation invariant by their definition; in dimension 



d = 3, also the tensors W{' , W 2 ' and W3' 1 are translation 
invariant due to the envelope theorems of Miiller [32]. Trans- 
lation invariance is important whenever a natural choice for 
the origin is not available. 

All MF and MT are homogeneous, i.e. they fulfill the 
homogeneity relation W^ S {\K) = X 3+r - v W^' s (K). Table I 
specifies the translation and homogeneity behavior of the MF 
and MT. 

Thus far, MT have been defined for (a) bodies with a smooth 
bounding surface (these may be non-convex) and (b) convex 
bodies which may have sharp corners and edges. The case 
of non-convex bodies with concave sharp corners or edges 
cannot be treated by the parallel body construction ("dilation") 



rl.l 



,W 3 



IM 3 = V with x{Br) = 1. 
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Homogeneity 
[unit] 


rank 


rank 1 


rank 2 translation 
behavior 


A 5 [m 5 ] 






Wq ' genuinely t. cov. 


A 4 [m 4 ] 


- 




W 1 ' genuinely t. cov. 


A 3 [m 3 ] 






W% '° genuinely t. cov. 
Wq Q t. invariant 


A 2 [m 2 ] 


W x 


- 


W^'° genuinely t. cov. 
W®' 2 t. invariant 
W\ Q t. invariant 


A 1 [m 1 ] 


W 2 




genuinely t. cov. 
W®' 2 t. invariant 
Wi Q t. invariant 


A" [1] 


w 3 




W3 Q t. invariant 



TABLE I 

Basic tensor valuations in 3D. The Minkowski functionals 

(SCALARS, RANK 0) ARE MOTION INVARIANT, THE MINKOWSKI VECTORS 
(RANK 1) ARE GENUINELY TRANSLATION COVARIANT. FOR RANK TWO, 

THE SPACE OF MINKOWSKI TENSORS DECOMPOSES INTO TWO 
COMPLEMENTARY SUBSPACES ACCORDING TO TRANSLATION BEHAVIOR 
(INDICATED BY THE LAST COLUMN): GENUINELY TRANSLATION 
COVARIANT AND TRANSLATION INVARIANT TENSORS . THE LATTER 
INCLUDE TENSORS OBTAINED BY MULTIPLYING THE SCALAR MF W v 
WITH THE UNIT TENSOR Q := e 2 + e 2 . + e| OF RANK TWO, WHERE 
ei, G2, e-i ARE VECTORS OF AN ORTHONORMAL BASIS OF R 3 . 



of the basic tensor valuations 2 [31, p. 150, line 23] Q P W^ 8 , 
that is 

<p{K) = VT P Q P W:- S {K) (12) 

with real coefficients ip]) s ' p that do not depend on the convex 
body K. The coefficients l p r j s ' p vanish unless r + s + 2p = 
2. Starting from eq. (12) and using the linear dependences 
among the basic tensor valuations, ip can be expressed in terms 
of linearly independent basic tensor valuations which form 
a basis of the corresponding vector space. The vector space 
of continuous, isometry covariant tensor valuations of rank 
two in R 3 has dimension 10. A particular basis of this vector 

2 2 2 

space consists of the six tensor valuations W ' , W, ' , W 2 ' , 

2 2 2 

Wg ' , W 1 ' and W 2 ' , which contain important independent 
shape information, and of the four tensor valuations QW V , 
v = 0, . . . , 3. A summary is provided in Table I. 

Since tp is continuous and additive on convex bodies, it can 
be extended as an additive functional to finite unions of convex 
bodies. For this additive extension, eq. (12) remains valid 
since the right-hand side is a linear combination of additive 
functionals. It should be emphasized, however, that although 
all these functionals are continuous on the space of convex 
bodies, they are not continuous on the space of finite unions 
of convex bodies, see the example in [6, Fig. 3]. 



without further assumptions (technically, such bodies do not 
represent sets of positive reach [21], [60]). An extension of the 
definition of MF and MT to finite unions of convex bodies is 
achieved by exploiting a property called additivity 

W^ S (KUK') = W^ s (K) + W^ s (K')-W^' s (KnK'), (9) 

if K, K' are all convex. In general, the union (K U K') of 
two arbitrary convex bodies K and K' is not convex while the 
intersection (K n K') is convex. Since W^ s are continuous 
functionals, eq. (9) can be used (see Groemer's extension 
theorem [21]) to define the MF and MT for bodies that are 
unions of a finite number of convex bodies (such sets are called 
polyconvex). 

MT of rank (r + s) with r + s > 1 are not completely 
linearly independent, i.e. they do not contain independent 
shape information [31], [61]. For rank-2 MT and d = 3 the 
linear dependences 

QW V {K) = vW^{K) + (3 - u)W^(K) (10) 

are valid for any polyconvex body K in M 3 and v = 0, 1, 2, 3; 
in particular, it follows that QW3 = 3W° ' 2 . Specifically, for 
v = one obtains a special case of the Gauss' theorem 

QWa = ZWl' 1 Wq = trWi' 1 . (11) 

These relations are special cases of [31, eq. (1.1)] or [61, 
eq. (1.5)]. 

Alesker's theorem [30] makes a strong statement about the 
completeness of the MT for the purpose of shape description. 
For the special case of tensors of rank two, it states that 
any isometry covariant, additive, continuous functional ip on 
general convex bodies in R 3 , taking values in the space of 
symmetric tensors of rank two over R 3 , is a linear combination 



A. MT of convex polytopes 

It is instructive to illustrate the MT for convex polytopes 
P and to point out similarities to the tensor of inertia. We 
use here the letter P, rather than K, to denote a body whose 
bounding surface is a polytope, i.e. consisting of piece-wise 
linear facets. For a polytope, the tensors W^'° characterize the 
distribution of mass if the cell is a solid cell (W n 2 '°), a hollow 

2 2 

cell (W-i 7 ), a wire frame (W 2 7 ) and a cell consisting of points 

2 

at the vertices only (W 3 ' ); in the last two cases, however, this 
distribution of mass is weighted with certain exterior angles 
(see Fig. 3). 

The tensor of inertia I, defined by 1^ = 
J K (— XiXj + Sij |x| 2 ) dV, is a measure of the mass 
distribution of a (homogeneous) body K, relevant for 
the relationship between a rotation and the resulting 
moment. As I is not translation-covariant, it is not a linear 
combination of the MT. However the simple relationship 
I(K) = -Wq'°(K) + tr (Wo'°(K^j Q holds for arbitrary 
K; as above Q is the unit tensor of rank two (metric tensor) 
of R 3 . This illustrates that W ' (K) is also a measure of the 
distribution of mass if K is a homogeneously filled solid. 
Similarly, W 1 ' (K) characterizes the mass distribution if 
K is hollow and bounded by an infinitesimally thin surface 
sheet. This interpretation of W^ and W 1 ' is valid for 
bodies K bounded by arbitrary surfaces (not just polytopes). 

For a polytope P, the tensor W2 ,0 {P) reduces to a line 
integral over the edges of the polytope (as the mean curvature 
vanishes on the flat facets, see also section II) and is hence 
related to a mass distribution if P is given by a wire frame with 
wires along the edges. However, imposed by the requirement 

2 An additive functional is called valuation. 
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(a) Wq'° - moment tensor solid (b) - moment tensor (a) (bl) (b2) 




wireframe vertices 



Fig. 4. Construction of a local parallel set. (a) Definition of the normal field 
over an arbitrary convex body K at its boundary dK. (bl-b2) Local parallel 
set Mi (P, rj) (illustrated for a polytope P) where only points are considered 
for which the normal direction is in a prescribed subset S of the unit sphere. 
(Here: t) = l d xS.) 



(e) W l : - normal distribution (f) W^ 2 - curvature distribution 




Fig. 3. Interpretations of MT for a polytope P: (a) Assuming homogeneous 

density, W 9 ' is the mass distribution tensor, (b-d) Contributions to surface- 

2 

integrated moment tensors W v ' with v = 1, . . . , 3 are concentrated on 
(3 — i/)-dimensional surfaces, (e-f) Moments of the normal distribution on 
dP. Contributions to normal moment tensors W^' 2 with v = 1, . . . , 2 are 
also concentrated on (3 — i/)-dimensional surfaces. 

of additivity, the weight of the wire cannot be uniform but must 
be proportional to the mean curvature along the edge (i.e. the 
dihedral angle). Similarly, the tensor W 3 ' (P) reduces to a 
sum of point contributions, as the Gaussian curvature G3 of 
P vanishes except at the vertices of the given polytope P. 
Due to additivity again, these vertices are weighted with the 
Gaussian curvature G3. 

B. Definition based on fundamental measure theory 

This section describes the alternative (and in some sense 
more fundamental) definition of MF and MT based on integral 
geometry and fundamental measure theory. The purpose of 
this section is to bridge the gap between the mathematics 
and physics literature on MF and MT. However, its content 
is not required for the numerical approaches to MF and MT 
described in section II. 

In integral geometry, the definition of MF and MT is based 
on so-called support measures (sometimes called generalized 
curvature measures) that can be thought of as local versions 
of the scalar MF [18], [20], [21], [62]. 

If support measures are available, then the MT for convex 
(or more general) bodies are obtained by integrating tensor 
functions with respect to these measures. Here we describe the 
approach for convex sets. The idea underlying the introduction 
of support measures for convex sets is to generalize the notion 
of a parallel set (or "dilation", see Fig. 4 for d = 2) of a convex 
body K in d-dimensional Euclidean space M d by a suitable 
local construction. 



A definition of the local parallel set that also applies 
to convex bodies K without smooth boundary is given in 
the following: We define Pa'(x) as the unique point in K 
which is nearest to a given point x G IR . This defines a 
continuous map Pa'(-) : R d — > K, x H- p/c(x). Then 
dfj-(x) := || x — p#T (x) || is the distance from x to K and 
njc(x) := (x — p/f (x))/dii'(x), for x G R d \ K, is an 
exterior unit normal of K at the boundary point pk (x) € dK 
(see Fig. 4). This construction achieves a definition of nearest 
points (reminiscent of the Euclidean distance map [63]) and 
surface normals that is also well-defined for points x whose 
nearest point on K is a sharp corner (where hence, the 
tangent plane is not well-defined and hence the conventional 
differential geometric definition of the surface normal does not 
apply). 

Now we shall derive a local Steiner formula and support 
measures, following ref. [21] to obtain local support measures. 
Then the definition of MF and MT follows directly as a special 
case [31], [61]. The intuitive idea underlying the definition of 
a local parallel set is described in two steps. First, we specify 
a region L C M d and some e > 0. Then we consider all points 
x G R d \ K which have distance dx (x) at most e from K 
and for which P.r:(x) G L. Second, if K has points (such 
as at sharp corners or edges) where the exterior surface unit 
normal vector is not unique, then it is natural to restrict the 
points in this local outer parallel set further by also requiring 
n/f(x) to lie in a prescribed set S C As an example, 

consider the polytope P in Fig. 4 (bl); for the definition of the 
local parallel set (shaded dark) it is necessary to specify which 
angular fraction of the wedges should be part of the local 
parallel set. This motivates the definition of the local parallel 
set by specification of the spatial region L and by the subset 
S C of normal directions, conveniently combined to 

T] = L x S C R d x where is the d-dimensional unit 
sphere. This two step procedure can be merged and slightly 
extended to the following general definition. 

For given e > and rj C M d x the local parallel set 

of K defined by 

M e {K,ri) :={xeM d \K :d K ( X ) < e, (p*-(x), n^x)) G rf\ 

(13) 

contains all points x G R d with < dfc(x) < e such that the 
pair (pa'(x), n^-(x)) G rj. The first condition restricts x to the 
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Fig. 5. (a) Illustration of the set of normal vectors on the surface of a 
polytope P. (b) Normal vector assigned to a facet F £ .^(P). (c) Line 
segment of normal vectors assigned to an edge F £ .Pi(P). (d) Spherical 
triangle/patch of normal vectors assigned to a vertex F £ J-"o(P). 



global outer parallel set (K \+)B e )\K (see Fig. 4 (bl)). The 
second condition restricts A4 e (K,r]) to those points x G R d 
of the global outer parallel set, where (p^(x), n^(x)) G rj. 

The volume of this local parallel set of a convex body K 
is Vd(Ai e (K, ij)). A fundamental result in integral geometry, 
known as the local Steiner formula [21], [62], states that the 
map e (->• Vd(A4 e (K, 77)) for all e > (V d is the d-dimensional 
Lebesgue-measure) is a polynomial of degree d, that is 

d-l 

V d (M € (K, 77)) = J2 e d - v K d - v A v (K, 77), (14) 

where n n :— n n / 2 /T(^ + l) is the volume of an n-dimensional 
unit ball and A V (K, rf), v = 0, ...,d— 1, are certain real 
coefficients that depend on K and 77, but not on e. 

For eq. (14) to be true for all e > 0, it is crucial that K is 
convex [64]. Equation (14) is easily confirmed (and evaluated) 
for a convex polytope P. In this case, the set (K WB e )\P can 
be decomposed in an elementary way into wedges over the 
faces F of P as indicated by Fig. 4 (bl). Let J- V (P) denote 
the ^-dimensional faces of the polytope P 3 and 



n(P,F) :={n P (x) G 



iti-l 



Pp(x) G relint F, x G R d \P} 

(15) 



the set of unit normal vectors assigned to F G F U (P) 
and relint F the relative interior of F 4 ; see Fig. 5. The 
contributions to V d (M e (K, 77)) coming from each of these 
wedges 

M e (P,F) :=M e (P, (relint Fx S^ 1 )) (16) 

3 T\ is the set of edges. In sec. II we use this notion for oriented edges, that 
is (non-oriented) edges are split into two oriented edges pointing in opposite 
directions. Here the set T\ contains non-oriented edges. Elsewhere it is stated 
explicitly. 

4 i.e. the interior of F w.r.t. the lowest-dimensional embedding affine hull 



can be calculated by a simple integration (known as Fubini's 
theorem or Cavalieri's principle [65], [66]). This shows that 
for 77 = L x S eq. (14) holds with 

1 



A u {P,LxS) = 



E 



FDL Jn(P,F)nS 

(17) 

where lo v = v k u is the surface measure of the (y — 1)- 
dimensional unit sphere and dTl 1 * is the Hausdorff measure 
of dimension v [67]. 

Since an arbitrary convex body K can be approximated 
by polytopes, eq. (14) can be derived by a continuity ar- 
gument. In fact, A v {K,rj) can be expressed as a linear 
combination of Vd(M e (K, 77)), for e M > pairwise different 
and /1 = 1 , . . . , d. One obtains an invertible d x c?-matrix 
equation (with the matrix entries e d ~ v Kd-i/) with v running 
from 0, ...,d — 1. Therefore, the properties of the local 
parallel volume Vd{M e (K, 77)) (in particular additivity, weak 
continuity) are also available for A v (K,rf) [62, p. 202]. In 
particular, K i-> A v (K, 77) is an additive functional for fixed 
77. That is, for convex bodies K and K 1 

A V {K U K', 7;) = A V (K, v) + A V {K', 77) - A V (K n K\ 77). 

(18) 

Furthermore, 77 H» (K, 77) is a non-negative measure for 
fixed K. The latter means that if 77^ C R rf x /j, G IN, is 

a sequence of mutually disjoint (measurable) sets, then 



E 



K{K, 77 M ). 



(19) 



A A K, [J 7 h 

This property is called er-additivity of A V (K, •). Weak conti- 
nuity of A„ means that for every sequence of convex bodies 

K(n) (a* S IN)) with K{u) ~^ K an( l ever Y continuous function 
/ : R d x ffi''- 1 -> [0, 00) the equation 



(20) 



lim / f(x,n)A l/ (K { ),d(x,n)) 

I— too J V 1 

= J /(x,n)A„(Jf,d(x,ii)) 



holds. Note that this does not imply A v (K( p ^rf) — s- A„(F, 77) 
as /i — > cxd for all measurable sets 77 C R d x 

In particular, A V {K,-) can be used to integrate functions 
over R d x $ d_1 . It is plausible that A U (K, •) is concentrated 
on the normal bundle N(K) := {(p^-(x), n^(x)) G dK x 
: x G H d \ K}. In other words, N{K) consists of all 
(x, n) G dK x such that n is an exterior unit normal 

vector of K at x. The measures A U (K, •), v = 0, . . . , d — 1, 
are called support measures and are determined as coefficient 
measures of the Steiner formula (14). They are local versions 
of the classical MF W v , since A V {K, R d x S^ 1 ) = V U (K) oc 
W d _ v {K). 

If K is sufficiently smooth, then 

A, (K, 77)- " ' 



i- f 1 {(x, njt(x)) G 77} G d _„(x) (U, 

" ^ (21) 
where G„(x) is the (1/ — l)-th (normalized) elementary sym- 
metric function of the principal curvatures of dK at x. That 
is in three dimensions, these are 1, the mean curvature and 
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Gaussian curvature, respectively. 1{-} is the characteristic 
function, which is evaluated to one if • is true and otherwise. 
For general dimensions and v < d — 1, eq. (21) holds for all 
convex bodies K. 

Having introduced the support measures as local versions 
of the scalar MF, it is easy to define the MT for a convex body 
K by 



1 Ud-v 

-!.«! 



x r n s A„ (K, d(x, n)) 



(22) 

hence we obtain ^ r v ' s {K) by integrating the tensorial function 
x r n s with respect to the measure A„(AT, •) over N(K) C 
R d x If if is a polytope, for d = 3 this yields equation 

(48) (up to a different normalization). If K is smooth, we 
obtain eq. (5) (up to a different normalization and indexing 
scheme), i.e. 

K' S (K) = . \ " ^ / x r n s G d _,(x) dA, 
r\s\oj d - u+s J dK 

^'V) = i / xW ' W 
r - Jk 

The notation $ M , r , s or for the MT in eq. (22) is 
preferred in the mathematical literature and differs from the 
notation Wl' s in eqs. (4-5) only by a different indexing scheme 
and a different normalization. In R 3 , i.e. for d — 3, the 
functionals 3>J; S and W£ ,s are related by 



W^(K)=C^/_ V (K), with 



C :-- 



for v = 1 , . . . , d, and 

W^(K)=r\^f(K). 



(24) 



(25) 



The additivity and continuity properties of the support 
measures immediately yield the corresponding properties of 
the MT. This approach also shows that if it is possible to define 
support measures for a class of sets, then the corresponding 
tensor valuations can be defined by eq. (22). 

Since the theory of support measures is well-developed [21], 
[62], the measure theoretic approach outlined above has some 
advantages over the differential geometric approach. 

As a simple illustration, let us explain why W\^ v is 
translation invariant for v = 0, . . . , d — 1. Observe that by 
translation covariance of the support measures 



/ xnA„ ( K&t, d(x,n) 

= [ (x + t)nA„ [K, d(x,n) ) 

jR d xS d - 1 \ / 

= / xnA, y [K, d(x,n) ] + 

JR d xS d - 1 \ J 

t / nA„ [K, d(x,n)J . 



(26) 



It is a basic property of the measures A U (K, •) that they are 
centered at the origin in the sense that J n A V (K, d(x, n)) = 
[32], which yields the assertion. 



A natural and useful extension that is suggested by general 
measure theory is to introduce local tensor valuations by 
restricting the integration on the right-hand side of eq. (22) 
to measurable subsets of R d x 

II. Bodies Bounded by Triangulated Surfaces 

We describe an exact algorithm for all independent scalar, 
vectorial and rank-2 MT of bodies bounded by piece-wise 
linear (i.e. triangulated) surfaces. Henceforth such bodies, 
convex or non-convex, are called polytopes P and their 
bounding surface is the triangulation T. Triangulations are 
commonly used as discrete approximations of smooth surfaces. 
The continuity of the MF and MT guarantee the convergence 
of the formula for the triangulations to the MF and MT of the 
smooth body. 

The formulae are derived for convex bodies with triangu- 
lated bounding surfaces by considering parallel bodies P e of 
convex polytopes P (that is P e has a continuous normal fields 
and finite curvatures for e > and a well defined limit as 
e — > 0). By application of the additivity relation these formulae 
are then shown to be valid also for bodies that are not convex. 
As the key results of this article — explicit formulae for the 
computation of MT of convex and non-convex polytopes — 
are summarized in table II. 

Consider a polytope P in 1R 3 with piece-wise linear bound- 
ing surface dP = J ' . Without loss of generality the linear 
facets may be assumed to be triangles. 5 The set of all 
triangular patches of dP is Ti, the set of oriented edges 
is T\ and the set of vertices is T§. We assume a doubly 
connected edge list (DCEL [68], also called half edge data 
structure [69]), that is, every edge which is shared between 
two triangles T and T" is a double-edge consisting of two 
oriented edges e (being part of T) and e' (part of T"), 
constituting an unambiguous assignment of each edge to a 
triangle. Each oriented edge is assigned to its previous edge 
e P revious and its next edge e ncxt . The remaining ambiguity in 
the edge orientation is lifted by requiring the triangle normals 
n T = (ep r0 viou S x e)/||e previous x e|| to point out of the body 
P. Thus we can uniquely assign to each oriented edge e a 
triangle T with vertices vi, v 2 v 3 (see Fig. 6) and its normal 
vector riT- 

The parallel body construction is illustrated by Fig. 2 (b). 
For an arbitrary body P the parallel body P e with thickness 
e > is defined as P e := P l±l B € . For a convex polytope 
P (whose bounding surface has a discontinuous normal field) 
the bounding surface dP e has a continuous normal field. The 
curvatures are patch-wise constant: G2 = G 3 = on the 
planar patches, G 2 = (2e) _1 and G 3 = on the cylindrical 
patches corresponding to polygon edges, and G 2 = 1/e and 
G3 = 1 /e 2 on the spherical patches corresponding to polytope 
vertices. For convex polytopes, the MT are defined as the 
surface integrals of eq. (5) evaluated on dP e in the limit e — > 0. 
The result thus obtained is equivalent with eq. (22), see also 
eq. (48). 

5 It is an important consequence of the additivity relation that the MT (in 
contrast to e.g. the texture tensor) do not change if flat polygonal facets 
are broken up into triangles. This is evidently also true for the algorithmic 
implementation described here. 
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Fig. 6. (a) Definition of geometric properties of a triangulated surface dP with doubly connected edge list (DCEL). Each edge e is oriented and uniquely 
assigned to a triangle T. The counter oriented edge to e is denoted e' and assigned to the adjacent triangle T'. An oriented edge e is unambiguously assigned to 
the previous edge e prcv ; ous and the next edge e ncx t. The normal vector ny is defined to point out of the body K, i.e. = (e prcv ; ous X e)/||e prcv ; ous X e||. 
The angle between two edges of triangle T at vertex v is denoted <j>^,. (b) Cross-sectional view along an oriented edge e. The normal vectors and n T i 
of the triangle T (that contains e) and T 1 span the angle a e £ (— n, n]. A concave edge has a negative angle a e . The figure also shows the definition of the 
local coordinate system used for the computation of W®' 2 . The basis vectors ri e , ri e and e are defined as e = e/|e|, ri e = (ny + njv)/|n;r + n T / 1 and 
ri e = h e X e. (c) Subdivision of a body P along a concave edge e to yield locally convex bodies. 







scalar measures 




Wo 






i E <C T ,n T )|T| 
T6J 2 


llj 


3 IdK dA 




| E |T| 

T6J 2 


w 2 






T5 E |e|Oe 


w 3 


U 9K G 3 dA 




3 E (*r- E 0?O 

v£/b TE7 2 (v) 






vectorial measures 






f K mdV 




E (y«(nr)*i see sec. II-D 

TGJF 2 




5 IdK Xi dA 




| E |T|(C T )< 




S IdK G2 X » dA 




f5 E |e|«o(C e )i 

e6Ji 




S IdK G 3 x » dA 




3 E (2i- - E ^)vi 

v£/b TE7 2 (v) 






tensorial measures 






Ik x » x i dv 




E {JT)ijk{n T )k, see sec. II-E 




3 IdK XiX J dA 




1 E (/t)« 




3 IdK G2 X * X J dA 




i E a e |e| ■ ((vf)ij + (viV 2 )y + (v|)<j) 
ee.Fi 




3 IdK G 3 X * X J d ^ 




| E ( 2n - E *?. J (v 2 ) y 




3 IdK n i n J dA 




| E |T|(n5,)y 

T6J 2 


(K 2 h 


3 IdK G 2 ninj dA 




^4 E l e (( a e + sina e )(n|)ij + (a e - sin Cf e )(h e )ij) 



Second column: MF and MT for bodies with smooth boundary dK. The mean and Gaussian curvatures on dK are G2 and G3, respectively. TTu'rd column: 
MF and MT for a triangulation. The set of facets of the triangulation T is T2, the set of oriented edges is J-i (in DCEL structure, see text) and the set 
of vertices To. The subset of triangles that contain the vertex v is denoted by J-^v). The nomenclature for triangulated surfaces is defined in Fig. 6. The 
vertices of an edge e or a triangle T are denoted vj, V2 and V3, respectively. |T| is the area of T 6 T2, Ct its center point (vi + V2 + vs)/3 and the 
tensors I? and are given in eqs. (33) and (35) and table III. C e is the center point of edge e; C e = (vi + V2)/2. k 6 {x, y, z}. (see also [58]). 



TABLE II 

MF AND MT IN 3D OF BODY K WITH SMOOTH BOUNDARY dK AND A BODY P BOUNDED BY A TRIANGULATED SURFACE dP. 
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Fig. 7. Subdivision of a non-convex body into convex sub-bodies, P = 
Pi U P2 U P3 U P4 U P5. Note that for the computation of MT the segments 
P2 and P4 need to be taken into account, though their volume measure is 0: 
Wo(P2) = W (P 4 ) = 0. 



A. Volume Wq 

The calculation of the volume of a polytope P can be 
transformed into a surface integral by Gauss' law (see eq. (11)) 
[70]. With divx = div(a;, y, z)* = 3 one obtains 

W (K) = [ dV = J I divxdF = \ f (x, dA) 

Jp 3 Jp 3 Jgp 

= tr - / xn dA = tr W^' 1 (27) 

JdP 

where dA = ndA denotes the oriented infinitesimal area 
element and {•, •) the scalar product. 

B. Surface area W\ and integral mean curvature W2 

The surface integral is a sum over triangles and is easily 
evaluated yielding the formulae in Tab. II. This result is 
independent whether P is convex or not. The surface area 
Wi(P) of dP is simply the sum of triangle areas. 

Expressing W2(P) as the limit of vanishing parallel dis- 
tance e of W 2 (P e ) of the parallel body P e , W 2 (P) = 
lim e ^o W^(P £ ), the contributions of facets vanish because the 
mean curvature of a flat face is zero. The contribution of the 
spherical patches S corresponding to vertices vanishes because 
the integral over a spherical patch S can be parametrized in 
spherical coordinates by ff s \ e 2 sin 9 dip d6 which vanish as 
e — » 0. The remaining contribution of the edges is located at 
the cylindrical patches of dP t and is given in polar coordinates 
by 

(28) 

where |e| is the length of edge e and a e the dihedral angle, 
i.e. the angle between the surface normals of the two facets 
adjacent to e. For a convex body, all edges have a dihedral 
angle < a e < tt; see also Fig. 6. Note that T\ is the set 
of oriented edges, i.e. the edge shared by two triangles is 
represented by two distinct oriented edges, which explains 
the factor 1/2 in front of the integral. 

Eq. (28) remains valid even if P is not convex, as is 
shown by exploiting additivity: A non-convex polytope P 
can be decomposed into a set of convex polytopes by cutting 




Fig. 8. Sketch of a vertex v with a spherical patch S of the parallel surface 
8P e . The interior angle in the triangle T adjacent to v is denoted tj>Z,. S is 
a spherical polygon. The jump angles coincide with the interior angles of the 
triangles. 

along the symmetric bisector planes of all concave edges 
(that is, — tt < a e < 0), see Fig. 7. For a concave edge e, 
the symmetric bisector plane is the plane that is spanned 
by e and the average of the facet normals of the two facets 
adjacent to e. By adding the contributions of all resulting 
convex bodies using the additivity relationship eq. (9), as 
outlined in Fig. 6 (c), one obtains the validity of eq. (28) 
for non-convex triangulated bodies. The sign of the dihedral 
angle a e £ (— 7r, tt] determines if the edge is convex (a e > 0) 
or concave (a e < 0). 



C. Integral Gaussian curvature W3 (Euler index x) 

As the point-wise Gaussian curvature G3 on cylinders 
and flat facets vanishes, only vertices of the triangulation 
(and their corresponding spherical patches on the parallel 
body) contribute to W3. For both a convex and a non-convex 
polytope P the point-wise Gaussian curvature G3 and the 
integrated Gaussian curvature W3 can be calculated by the 
well-known simple sum over angle deficits at surface vertices 
in eq. (30), derived below, and also given in [71], [72]. The 
non-convex case is treated by exploiting additivity. 

The Gaussian curvature contribution of the vertices v£ Jj 
is derived by the Gauss-Bonnet-formula 

f G 3 dA = 2n- 4>t- I fc s ds > ( 29 ) 

where ^(v) is the subset of triangles adjacent to vertex v 
and S denotes the spherical patch on the parallel surface dP t . 
For all e > each spherical cap S C dP e can be uniquely 
assigned to one vertex v; dS its oriented boundary curve and 
kg the geodesic curvature along dS. At the corners of dS, the 
discontinuity of the tangent vectors is characterized by jump 
angles <f>X, (see Fig. 8) which for all e > coincide with the 
interior angles of the triangle T at v [73], see Fig. 8. The 
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X, z 
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y, z 


(xyz,0,0f 


X 



TABLE III 

Auxiliary functions used to convert the volume integrals of 
and W 2 ' 



Fig. 9. The Gaussian curvature G3 at a saddle vertex is obtained by a virtual 
decomposition of P at v into three polytopes with vertices v M , fj, = 1,2, 3 
using the additivity property of MT. 



geodesic curvature k g vanishes almost everywhere along dS, 
because dS are great circle arcs on the spherical patch and 
the adjacent cylindrical patch and are thus geodesies; hence 
the integral J gs k g ds vanishes. 

As a consequence, J s G3 cL4 is constant for all e > 0. 
Equation (29) therefore yields a definition and an explicit 
formula for W 3 (P) as a sum of the local contribution wz(P, v) 
at a vertex v 



can be explicitly written as 



w 3 (P) = ^ ]T W3 (p,v) 



3 ^ 



2tt 




At a concave vertex v, a polytope P can always be decom- 
posed into three separate bodies (one of vanishing volume) 
that have convex vertices in lieu of v. It is easy to validate 
eq. (30) at concave vertices by using the additivity relation in 
eq. (9) (see Fig. 9). 



D. Center of mass W ' and curvature centroids W^'° 

The Minkowski vector W ' corresponds to the center 
of mass of P multiplied with its volume (assuming P is 
homogeneously filled with material of constant density.) The 
components of this vector may be computed by transforming 
the volume integral into a surface integral using Gauss' 
theorem 



(7 T )ife(nT)fc 



(32) 



with k as listed in table III. (k is not a summation index.) 
\T\ is the surface area of T. The It in eq. (32) are tensorial 
integrals over the individually parametrized triangles with the 
three vertices v M , fi = 1,2,3 

1 l-a 

It = 2\T\ J da J db [ Vj + a(v 2 - v x ) + b(v 3 - Vl )] 2 . 


(33) 

The components of the auxiliary functions (fi)k are selected 
entries of the tensor It or zero. It can be written in terms of 
the triangle vertices v M and triangle centers Cj of T as 



(34) 



The remaining integrals W„'° with v = 1,2,3 are evaluated 
similarly to the integrals W v W%'° (see below). The integrals 
W®' 1 (v = 1,2,3) involving surface normals vanish for 
arbitrary bodies (with closed bounding surfaces). 



E. Volume integral W 

The volume integral Wq ,0 (P) can be computed in a similar 
way as W 1 '°(P). Using 

1 l-a 

J T = 2\T\J da J db [vi + a(v 2 - Vi) + 6(v 3 - v x )] S 







(w^°(P)J = y x, dV = J (fi, dA) (31) the components ij of the tensor may be expressed as 

p 9P / ., r, \ » ^ . . . 



(35) 



with functions that satisfy divf; = x^. The vector- valued 
auxiliary function can be chosen for each index i indepen- 
dently and the index i denotes the index of Wr}' , which is 
evaluated. For the particular choice of ii given in table III, this 



(WS'°{K)) = (Jr)«*(nr)*. (36) 



TGP 2 



Again, the index k is not a summation index but rather the 
index specified in table III. This derivation applies equally to 
convex and non-convex polytopes P. 
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F. Surface integrals W®' 2 and W 2 ' 2 

2 

The computation of W 1 ' results in a simple sum of 
integrals over triangular facets, resulting in the formulae in 
Tab. II, both for convex and non-convex bodies. 

The tensor W 2 ' 2 is calculated by a parallel body construc- 
tion, first demonstrated for convex bodies. Consider a convex 
poly tope P, and the corresponding parallel body P e . The 
integral over the parallel surface is split up into integrals over 
flat facets, cylindrical patches and spherical patches. Out of 
these only the cylindrical edge segments contribute, for the 
same reasons as for the scalar measure W 2 . The remaining 
contribution is calculated for e — > using the following rep- 
resentation for the normal vectors on the cylindrical patches. 
Given an edge e with facet normals ny and tiy of the adjacent 
triangles. One obtains (also representing a special case of 
eq. (48)) 



W'r 



0.2 



(K) = 1 £ |e| 

ee.Fi 



b/2 



/ 

-a e /2 



n 2 M. 



(37) 



To compute the integral on the right hand side we define the 
orthogonal unit vectors e = e/|e|, ri = (n eT + n eT ,)/|n eT + 
n er , | and ri = e x ri. For a given edge, n(i?) can be written 
as n = cos $h + sin ■& ri. In this basis, the integral over 
n 2 evaluates to (1/2) ((a e + sina e )ri 2 + (a e — sina e )ri 2 ), 
see Fig. 6. This yields the formula in Tab. II. The validity 
of this formula for non-convex bodies follows from the same 
additivity arguments as for W 2 . 



2 2 

G. Curvature-weighted surface integrals W 2 ' and W 3 ' 
The mean and Gaussian curvature weighted surface integrals 



W 2 ' a and Wg ,u over position vectors can be evaluated as the 
limit e — > of the parallel body construction, for convex 
bodies. The validity for non-convex shapes follows from the 
analogous construction as for W 2 and W3 (see tab. II). 



r 2,Q 



H. Open bodies, labeled domains and Minkowski maps 

The analysis presented so far has been derived for compact 
bodies in R 3 with a closed bounding surface - and inherits 
strong robustness from its integral nature. For some analyses, 
the requirement of closed bodies is too stringent. For example, 
experimental datasets of percolating or periodic structures, 
both of which extend infinitely through space, always represent 
finite subsets of the structure with components that traverse 
the dataset boundaries. Similarly, an analysis of a periodic 
model may be restricted to a translational unit cell, see Fig. 10. 
Furthermore, a local MT analysis, termed a Minkowski map 
[6], [13], can be useful to quantify variations throughout the 
sample. For a Minkowski map, a grid is superposed on the 
body K, and the MT are computed separately for each grid 
domain L. Such Minkowski maps can be useful to analyze 
spatial heterogeneity of anisotropy or orientation at the length 
scale given by the size of L. In these situations, the MT are 
computed for the subset L(ldK of the whole bounding surface 



that is contained in a box L. In general, LC\dK is not a closed 
surface even if dK is. 

It is evidently possible to take the subset K n L of K 
and consider d(K n L) as the bounding surface. However 
this introduces bounding surface patches (e.g. solid/void 
interfaces if K is a porous medium) that are not part of the 
bounding surface dK of K. For physical analyses one may 
want to avoid such boundary effects, i.e. not consider the 
contributions of these additional bounding surface patches. 
This motivates the introduction of MF and MT for open 
bodies, i.e. bodies without a closed bounding surface (see 
Fig. 10). 

In lieu of an attempt to define MF and MT for open bodies, 
we define a domain-wise analysis of MF and MT. Consider 
a decomposition of the surface dP of a triangulated body P 
into m domains, or patches, D a (with a = l,...,m) such 
that 



dp = |J d„ 



(38) 



and consider these domains labeled by labels a = 1, . . . , m. 
Triangles are uniquely assigned to a label, but edges and 
vertices of the triangulation can be shared between several 
domains. Specifically, the domains D a could represent a 
decomposition of P into patches contained within the grid 
domains of a three-dimensional lattice (See Fig. 11). 

Contributions of the facets can be uniquely assigned to D a , 
but the contribution of edges and vertices on the boundaries 
of a patch D a needs to be divided between a and the label of 
the adjacent domain a' (See Fig. 11). 

For W 2 S , the contribution of the dihedral angle at an edge is 
equally divided between the labels of the two adjacent triangles 
(Note that this is naturally taken account of by the use of 
oriented edges in the doubly-connected edge list, discussed 
above). For W^' s the division of the contribution of the interior 
vertex angles to the integral Gaussian curvature measures W^' s 
is less straight-forward. An intuitive way, that is also consistent 
with global integration over all labeled domains, is provided 
by the use of label factors. The label factor fr> a (v) of domain 
D a at vertex v is defined as 



E Pt 

TGF 2 (v) 



(39) 



where J- 2 {D a ) is the set of all triangles labeled with label 
a. Hence fo a (v) is the sum of angles at v of those triangles 
adjacent to v and are labeled a divided by sum of these angles 
of all adjacent triangles. It is easy to see, that 



W 3(v) = £/A,(v)W3(v), 



(40) 



a=0 



for a vertex with m adjacent labels and W 3 (P) = 

EvGFo W 3(v). 

For the volume tensor Wq'° a label-wise analysis is only 
well-defined if the body K is subdivided (and not only the 
bounding surface dK). 
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Fig. 10. (a) A compact (but not convex) body K that corresponds to a 
translational unit cell of a periodic body and a triangulation of it. Its bounding 
surface consists in the translational unit as well as the flat "end caps" that seal 
it. (b) A surface portion representing a translational unit and its triangulation. 
The body K is the volume to one side of the surface and forms a connected 
periodic body. The surface and its triangulation extend beyond the translational 
unit cell indicating that the surface (and the triangulation) is periodic. 



(b) 





Fig. 11. (a) A triangulated surface may be decomposed into several domains 
by assigning a domain label to each triangle; also open boundaries are 
possible, i.e. triangle edges without adjacent triangles (thick dashed line), 
(b) Label factors at vertices which are adjacent to triangles of more than one 
domain. The Gaussian curvature contribution is weighted with the interior 
angles belonging to each domain. 



where the pi(n') = f gK S(n — n')dA That is, p\ is area- 
weighted density of normal vectors. It is easy to see that an 
uniform distribution on $ 2 is equivalent to an isotropic tensor. 

For example, if K is a sphere, then p\(n) is constant and 
/3°' 2 = 1 as expected. For the rectangular box [0, a x ] x [0, a y ] x 
[0,a 2 ], the function pi(n) is concentrated at delta peaks, 

5(ei - n) + 5(e. i + n) 



Pi( n ) 



i=x,y,z 



(43) 



The resulting anisotropy measure is a z /a x for a x > a y > a z . 
It is instructive to express the second translation-invariant 

2 

MT W 2 ' by a distribution of normals and curvatures. The 
density 



P2( 



n', G' 2 ) = [ S(n- n')S(G 2 - G' 2 )dA (44) 

JdK 

gives the sum of the area of all surface patches that have 
normal direction n' and mean curvature GL 



W°> 2 (K) 



G 2 / p 2 (n,G 2 )n©ndf>dG 2 . (45) 



s 2 



If the function p 2 can be written as a product 
p 2 (n, G 2 ) — p 2 (G 2 ) pi(n), the anisotropy characteristics /3?' 2 

2 

and (i 2 defined as the ratio of the smallest to the largest 
eigenvalue of W 2 ' are identical. In this sense, f3^ 2 provides 
a higher order anisotropy measure that quantifies anisotropy 
of the curvature distribution. 



III. Anisotropy Measures 

Based on MT, robust measures of anisotropy can be defined 
that are sufficiently sensitive to capture subtle anisotropy ef- 
fects and that are applicable to a wide range of microstructures. 
The usefulness and versatility of this approach is demonstrated 
by two examples representing different types of structures - a 
cellular partition and a network structure. 

A rank-2 tensor is defined to be isotropic if and only if 
it is proportional to the unit tensor Q, i.e. its eigenvalues 
are all equal. Deviations from isotropy are measured by the 
anisotropy index /3£' s , which is the ratio of extremal eigenval- 
ues of the tensor W^ s . For example, let £ M (|£i| < £ 2 | < I63D 



be the eigenvalues of W l ' 



jQ.2 



then the anisotropy index is 
$1 



6 



G [0,1]. 



(41) 



By definition, this quantity is dimensionless, continuous and 
rotation invariant. The value of 1 indicates perfect isotropy, 
and smaller values indicate anisotropy. For anisotropic bodies, 
it is sometimes also useful to consider 7°' 2 = |£ 2 /£3|. 

These quantities can be easily interpreted for the translation 
invariant tensors W^ 2 and W 2 ' 2 . We can write W^' 2 equiva- 
lently to eq. (22) as the second moment of the distribution of 
normal vectors (with the density pi{n)) on the unit sphere S 2 
as 



pi(n) n n dfl, 



(42) 



A. Alignment of Actin Biopolymer networks under shear 

Biopolymer networks made of actin or collagen fibers are 
important structural elements in biological tissue that act as a 
scaffolds and provide stiffness and mechanical stability [75]- 
[77]. Of current interest is the relationship between fiber 
arrangement and alignment on the one hand side and elastic 
or visco-elastic properties on the other. This relationship can 
be probed by shear-experiments with confocal microscopy 
providing real-space structural data [74]. We now demonstrate 
that the degree of alignment and of structural anisotropy of the 
fiber network is well-captured by a MT analysis. 

The data sets analyzed here represent actin fiber networks 
reconstituted from rabbit actin biopolymer networks with actin 
concentration of 1.2mg/ml cross-linked with filamin A. These 
are imaged using confocal microscopy for different shear 
deformations, see the explanation in Fig. 12. The data sets 
are the same as those analyzed in [74]. The gray-scale data 
set is converted into a binary data set with 1 corresponding to 
actin and corresponding to the surrounding fluid by standard 
threshold segmentation with threshold I c 6 . The medial axis 
of the 1 phase is computed using distance-ordered homotopic 

6 The threshold I c is chosen such that only the brightest and hence thickest 
fibers are retained. For a given segmentation threshold I c the integrated 
intensity of all voxels of the fluid phase is A = (J^ I{p))~ ^2* I{v) where 
I(p) is the intensity of voxel p in the original intensity data set, ^2 me sum 
over all voxels of the data set and J^* the sum over all voxels of the fluid 
phase, i.e. those voxels that are set to by the segmentation process. The 
values of I c chosen here correspond to 0.95 < A < 0.99. 
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Fig. 12. Anisotropy measures l//3^' 2 and lfy®' 2 of actin network as 
function of shear. The ratio of largest to the smallest eigenvalue 1 / /3®' 2 grows 
with increasing shear indicating that the fibers become increasingly aligned. 
The ratio of the largest to the intermediate eigenvalue ' 2 remains close 
to unity as expected since the fibers are essentially one-dimensional lines, 
inflated to approximately cylindrical tubes. The insert shows the alignment 
angle <f> of the eigenvector corresponding to the minimal eigenvalue £i of 
W-^ ' and the direction of applied shear; its decay to indicates that the 
network aligns with the shear direction when large shear is applied. The 
error bars in both plots are the standard deviation of the distribution of 
the quantities when analyzed for different segmentation parameter between 
0.95 < A < 0.99. The illustrations above sketch the experimental setup 
with confocal microscopy images corresponding to a fiber network at shear 
e := Al/h = and e = 2 (the small confocal microscopy images are taken 
from [74]). 



thinning [78], [79] and is used as the one-voxel thick line 
representation of the actin fiber network. Conversion to a 
triangulated representation is obtained by using the Marching 
Cubes algorithm [80]. For more details of the analysis of 
biopolymers see ref. [81]. 

Typically only a subset, or observation window, of the 
structure is available for analysis. Therefore we assume that 
the network is homogeneous and a sufficiently large but 
finite subset is accessible which is assumed to represent the 
entire sample. The derived measures /3^' 2 quantify the intrinsic 
anisotropy, i.e. their values do not depend on the size, aspect 
ratio or position of the observation window (for sufficiently 
large windows). 

Figure 12 shows l//?"' 2 and I/7J' 2 evaluated on the whole 
network (that consists essentially in a single component; only 
the outer boundary layers of the confocal data are clipped). 
It shows that the distribution of normal directions of the 
fiber bounding surface becomes less isotropic with increasing 
shear, indicating alignment of the fibers. The angle between 
the eigenvector to the minimal eigenvalue £1 (corresponding 
approximately to the dominant tangent direction) and the 
direction of shear decreases to 0, indicating the alignment of 
the fibers with the direction of shear. This is commensurate 
with the results published in [74], that extracted a distribution 
of tangent directions and used these to quantify alignment. 



The eigenvalue ratios of the translation-covariant tensors 
W 2,0 (and of the tensor of inertia) capture different aspects of 
the anisotropy of a shape compared to the translation invariant 
tensors W®' 2 , see also section I-A. However, usefulness of 
the translation-covariant tensors depends on whether or not a 
natural definition of the origin is available for the system. For 
example, for the analysis of cellular shapes one may choose 
the center of mass W ' /Wq or the corresponding curvature 
centroid W^'°/W v as the origin. Especially for percolating or 
periodic bodies, for which the analysis is always restricted 
to a finite window of observation, the choice of origin is 
often not naturally determined. An additional problem for such 
structures is that the measures /3 2 ' derived from translation 
covariant tensors W 2,0 , as opposed to the translation-invariant 
measures /3°' 2 , crucially depend on the shape and size of the 
window of observation. 

The analysis of alignment of biopolymer networks illus- 
trates the potential of the MT approach for structure character- 
ization of cellular and porous materials, and demonstrates its 
applicability to voxelized experimental data. The MT approach 
can shed light on systems with a similar spatial structure that 
exhibit subtle anisotropy effects, such as fibrous biological 
materials [58], porous materials [4], and metal foams [82]. 

B. Anisotropy of free-volume cells of random bead Packs 

Granular media are a physical system where geometry plays 
an essential role in determining the physical properties, such 
as flow or packing properties. The geometric structure to be 
characterized is substantially different from the above example 
and consists in an ensemble of grains. These grains do not 
fill space without gaps, and each grain is associated with its 
surrounding region of space, often through a construction of 
the Voronoi diagram. 

The distributions of volumes of the Voronoi cells of dis- 
ordered jammed sphere packs with packing fractions from 
0.55 (random loose packing) to around 0.64 (random close 
packing, RCP) have attracted interest for the study of granular 
systems [83], motivated by a possible statistical mechanics 
description of granular systems [84], [85]. Aste et al. [83] 
used the volume distribution of Voronoi cells to estimate 
configurational entropy in static packings. 

Experimentally, RCP is the maximal packing fraction ob- 
served in disordered monodisperse sphere packings. However, 
for non-spherical particles this limiting packing fraction has 
a higher value: anisotropic bodies such as M&M® choco- 
late candies pack more tightly in the disordered phase than 
spherical particles, reaching packing fractions of 4> « 0.71 
for spheroids and <fi rs 0.735 for general ellipsoids [86], [87]. 
Higher order shape measures of the Voronoi cells, such as their 
anisotropy, have only recently been investigated, demonstrat- 
ing characteristic anisotropy distributions in Voronoi diagrams 
of monodisperse spheres [1], It has been conjectured that 
ellipsoids are, due to their anisotropic shape, making more 
efficient use of the anisotropic accessible volume in a static 
packing and hence pack denser than spheres [1]. 

Experimental and simulated random close-packed configu- 
rations of approximately 10 5 — 10 6 monodisperse spheres were 
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Fig. 13. (a) Voronoi diagram of jammed monodisperse sphere pack, (b) The 
same Voronoi diagram, but spheres are replaced by ellipsoids with half-axes 
o-l = a,2 < oz, aligned along the eigendirections of Wq°. Colors represent 
the ratio of the shortest and longest axis of the ellipsoid. The relation of 
ai/a-j to /3 '° is given in Fig. 14. 



generated using three different protocols [88]-[90]. The center 
coordinates of all beads (in the experimental data sets) were 
extracted by X-ray microtomography and a deconvolution 
technique [91], [92]. The Voronoi cells of the sphere centers 
are computed using the program qhull [93]. For a point p of 
a set V of points, its Voronoi cell is the convex polytope that 
contains all points of R 3 closer to p than to any other point 
in V [94]. 

Figure 13 shows a subset of a experimental bead sample 
with packing fraction 0.58 on the left panel. The wireframe 
illustrates the Voronoi diagram of these spheres. On the right 
hand side, spheres were replaced by ellipsoids with the same 
anisotropy measure /3q'° and the half-axes a\ — a2 < 03. CI3 
is aligned along the eigendirection of the largest eigenvalue 

2 

of W Q ' of the corresponding Voronoi cell. The use of MT, 
and their eigenvalues and eigendirections, hence provides an 
efficient way of "fitting" ellipsoids to given convex cells. 

This substitution of spheres by ellipsoidal particles (allow- 
ing slight overlap) leads to packings with (j> ss 0.72, consistent 
with experimental and numerical observations for these kinds 
of packings [87]. 

This example has demonstrated the ability of MT to quan- 
titatively characterize the subtle, but important, anisotropy ef- 
fects in packings of granular particles. A quantitative analysis 
of this system, and of the relevance of anisotropy for an 
interpretation of maximally obtainable packing fractions in 
random granular media, is presented in ref. [1]. 

Voronoi volumes in granular matter can be connected to the 
configurational entropy [95]; the MT analysis goes beyond 
the description of volume alone. MT can be used to test 
local or global isotropy assumptions [1]. This is important 
for systems with inherent anisotropy, such as dense granular 
flows under gravity [96] or shear-zones in granular matter [97]. 
Furthermore, MT (of rank 4) have been proven useful for the 
precise determination of a crystallization onset in granular 
systems [7]. This demonstrates the relevance of advanced 
shape characterization for the study of granular matter. 

The characterization of cellular shapes is also relevant 
in various other contexts, e.g. for the relationship between 
structure and dynamics in glass-forming liquids [98] or for 
packing entropy of the hard micellar cores in supramolecular 
micellar materials [99], and for the understanding of physical 
properties of foams, both dynamic [100] and static [41], [101]. 



IV. Conclusion 

In this article, we have described linear-time algorithms to 
compute Minkowski tensors of three-dimensional triangulated 
bodies. We have shown that robust measures of intrinsic 
anisotropy can be derived from these tensors. These have 
been used to quantify the anisotropy of an experimental model 
system for granular matter and of confocal microscopy images 
of sheared biopolymer networks. 

The availability of this algorithm is a prerequisite for 
applying MT analysis more widely for the quantitative de- 
scription of anisotropic structural features in physical systems. 
Further applications of this algorithms in studies concerning 
anisotropy are discussed in refs. [2], [4] ranging from the 
analysis of Turing patterns, foams and cellular systems to 
defect detection in crystals. Rank-4 tensors have already been 
used to detect crystallization in jammed packings [7]. 

More and more 3D real space become available in many 
domains of science. This paper provides a efficient method 
to extract basic morphological information from experimental 
3D data. For the design of optimal spatial structures (see for 
example refs. [102], [103]) the understanding of the interplay 
between morphology and physical properties is essential. 

Software 

Software to compute the Minkowski tensors is freely avail- 
able at our website, 

www.theoriel.physik.fau.de/karambola 
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Appendix 

A. Specific examples 

For some simple shapes the rank-2 MT can be calcu- 
lated analytically by using explicit surface parametrizations 
and expressions for surface normals and principal curvatures. 
Specifically for a sphere of radius R centered at the origin, 
one obtains 



(46) 



and, for v = 1, 2, 3 and r + s = 2, 



4-7T 4-7T 

W v = —R 3 -", W r v s = —R 3 ~" +r Q (47) 



with the unit tensor Q. 
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For a convex polytope P, we write for the set of 

/j-dimensional faces of P, \i = 0,1,2, that is, J'o(P) is the 
set of vertices, F\(P) is the set of (non-oriented) edges, and 
T 2 (P) is the set of faces. If F € F^(P), then we denote 
by n(P, F) the set of exterior unit normal vectors of P at 
F, which is a (2 - fi) -dimensional subset of the unit sphere 
$ 2 . Then we obtain, as a special case of general formulas in 
section I-B, 

W:- S (P) = l E / x r W 3 -"(dx) / n"H v -\dn) 

,JF Jn(P,F) 

(48) 

with v = 1,2, 3. For the notation see the main text at eq. 17. 

For a cuboidal box of dimensions a x x a y x a z aligned 
with the coordinate axes and centered at the origin eq. (48) 
yields W = a x a y a z , W\ = ^{a x a y + a y a z + a z a x ), W 2 = 
f (a x + a y + a z ), W 3 = ^f, W^'° = and all MT of rank 
two are diagonal matrices with the following entries 

1 

:di djCik, 
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G 



(W 2 % 



— (af + 3a 2 ( aj + a fc )) 



where {i, j, fc} = {x, y, z} and permutations thereof. 

A torus centered at the origin with major radius Ri and 
minor radius R 2 < R\ can be parametrized by x(u, v) = 
{cos(u) (i?i + i?2 cos(w)) , sin(w) (i?i + i? 2 cos(u)) , R 2 sin(w)} 
with a, (3 G [0, 2tt). The scalar functionals are explicitly 
given by W = 2tt 2 R 1 R 2 2 , W x = ^R\Ri, W 2 = ^Ri 
and W3 = and the vectorial measures are W^'° = 0. The 
tensors of rank two are diagonal with degenerate eigenvalues 
{Wl' s ) xx = (W^ 8 ) yy ; the entries are given by 



Fig. 14. Eigenvalue ratio of the smallest and largest eigenvalues £ m i n and 
§max of the MT Wl' s of an ellipsoid with radii l x = 1 and l x = 1 > l y > l z 
as function of r = l z /l x - Each symbol in the main plot represents data (hardly 
distinguishable) for three different intermediate radii l y = 0.1,0.5,0.9 
indicating that for these four tensors the minimal to maximal eigenvalue 
ratio is approximately the same for all values of the intermediate radius. The 
solid curves are fits to the data giving /3o'° ~ 1.210r 3 — 0.235r 2 + 0.024, 
= r 2 , pf'° Ri /3°' 2 Ri -0.366r 3 + 1.222r 2 + 0.139r. The insert shows 
the eigenvalue ratio of the tensor W 2 ' as function of l z /l x . In contrast to the 
above four tensors, this ratio depends strongly on the value of the intermediate 
radius l y . In particular, for l z = the eigenvalue ratio only becomes zero if 



the intermediate radius is also l y 
ratio converges to 0.5 for l z /l x - 



0. For the maximal l y 
> 0. 



1 the eigenvalue 



minimal to maximal eigenvalue ratio of the MT of rank two 



of ellipsoids with 
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the surface integrals all result in elliptic integrals and cannot 
be expressed in closed form. However, the scalar MF Wo is 
Wo = ^-l x lyl z - The MT W 2 ' is diagonal with 

(W 2 '°) u - ^l%l k (49) 

where {i,j,k} is {x, y, z} and permutations thereof. The 
integration of all other tensors is easily numerically ob- 
tained by use of the ellipsoid parametrization x(u,v) = 
{l x cos(u) sin(w), l y sin(it) sin(7j), l z cos(w)} which yields ex- 
plicit expressions for the metric tensor of the ellipsoidal sur- 
face, the normal vector, and the mean and Gaussian curvatures. 
These are readily integrated numerically. Fig. 14 shows the 
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